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Wettability stabilizes fluid invasion into porous media via nonlocal, cooperative pore Ailing 
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We study the impact of the wetting properties on the immiscible displacement of a viscous fluid in disordered 
porous media. We present a novel pore-scale model that captures wettability and dynamic effects, including the 
spatiotemporal nonlocality associated with interface readjustments. Our simulations show that increasing the 
wettability of the invading fluid (the contact angle) promotes cooperative pore Ailing that stabilizes the invasion, 
and that this effect is suppressed as the flow rate increases, due to viscous instabilities. We use scaling analysis 
to derive two dimensionless numbers that predict the mode of displacement. By elucidating the underlying 
mechanisms, we explain classical yet intriguing experimental observations. These insights could be used to 
improve technologies such as hydraulic fracturing, CO 2 geo-sequestration, and microfluidics. 


PACS numbers: 47.54.-r, 47.56.+r, 47.20.-k, 47.55.-t 

Fluid-fluid displacement in porous media is important in nat¬ 
ural and industrial processes at various scales, from enhanced 
energy recovery, CO 2 geo-sequestration, groundwater contam¬ 
ination and soil wetting and drying, to dyeing of paper or tex¬ 
tiles and microfluidics. Fluid displacement is governed by the 
interplay between quenched disorder, short-range cooperative 
effects and long-range pressure screening, which depends on 
a large number of parameters, including the wettability—the 
relative affinity of the fluids to the solid. Consequently, the 
displacement patterns can range from a stable, compact front 
to highly ramifled with preferential flow paths (Angers) d. 
Fluid invasion is a member of a broad class of problems charac¬ 
terized by competitive domain growth and nonlinear interface 
dynamics, including magnetic domains, biological Aims and 
flame front propagation d. The interface evolution in these 
systems is often modeled as a competition between the energy 
associated with the interaction between phases and constraints 
arising from disorder; the relative importance of the two can be 
tuned by properties such as wettability in fluid displacement or 
local random interaction flelds in magnetic domains 0. Un- 
derstanding the impact of wettability on fluid invasion—the 
topic of this Fetter—is therefore relevant to a wide range of 
phenomena of scientiflc and technological importance. 

Immiscible displacement can be classifled according to the 
wettability into drainage or nonwetting invasion, where the dis¬ 
placed fluid preferentially wets the solid (contact angle 0 < 
90°, measured through the defending fluid), or imbibition of 
a wetting fluid {0 > 90°). Intensive research has provided 
basic understanding of drainage, identifying different invasion 
behaviors and explaining their dependence on the flow veloc¬ 
ity, fluid viscosities, interfacial tension, and the degree of pore- 
scale disorder ((SO and the references therein). Increasing 0 
was found to stabilize the displacement and reduce trapping in 
forced and gravity-driven drainage experiments Ea. 

In contrast, relatively few works have studied imbibition, 
mostly for the stable case of a more viscous invading fluid (H 
Biini. For unstable viscosity ratios, experiments showed a 
marked difference between viscous Angering in drainage and 
more stable patterns with thicker Angers in imbibition ifTTI . 
Stabilization was also captured in simulations which intro¬ 


duced viscous effects stochastically ifT^ and in lattice Boltz¬ 
mann simulations ca. Quasi-static simulations (neglecting 
dynamic effects) illustrated that increasing 0 enhanced the oc¬ 
currence of a nonlocal, cooperative pore Ailing mechanism, 
resulting in a compact pattern ifffil . These intriguing results 
were only recently explored systematically by experiments in 
which the wettability was altered while keeping the same fluid 
pair ca. The authors demonstrated that increasing 0 stabilized 
the displacement, leading to a compact front in slow imbibi¬ 
tion despite the high, unfavorable viscosity ratio ina. Many 
of these important observations remain unexplained, primar¬ 
ily because of nonlocal pore Ailing dynamics, that is inaccessi¬ 
ble experimentally and not well-characterized by existing mod¬ 
els El da. In this Fetter, we present a novel pore-scale model 
that exposes the competing effects of wettability and flow rate, 
thereby explaining the aforementioned observations. 

We develop a two-dimensional, discrete model of immis¬ 
cible displacement in a random medium with fluids of arbi¬ 
trary viscosities and contact angle. Our model is briefly de¬ 
scribed below, and in further details as Supplemental Mate¬ 
rial El- A mechanistic description of the displacement dy¬ 
namics with both capillary and viscous forces is obtained by 
combining two modeling approaches: (a) grain-based dm, 
resolving meniscus stability from pore geometry; and (b) pore- 
based (21 [13, resolving fluid pressures and fluxes from the 
pore topology and geometry. Through consideration of vis¬ 
cous dissipation, our model captures the nonlocal nature of 
interface dynamics: the effect of local pore invasion on the 
interface conflguration elsewhere, the disparate timescales of 
pore Ailing and bulk flow ||20| 1^ . and the associated mech¬ 
anisms of pressure screening (221 EH and interface readjust¬ 
ments (20l[24l. These mechanisms are crucial even in slowly- 
driven systems, limiting the volume invaded in a single event 
(avalanche) EOlEll and increasing trapping Ea, and, for 
the conditions considered here—high disorder and porosity, 
small throat to pore size ratio, and large viscosity ratio— 
have a stronger impact on the displacement than other mech¬ 
anisms such as film flow, snapoff, and contact angle varia¬ 
tions iniEiiia. Furthermore, modeling contact line and contact 
angle dynamics is strongly debated, and requires consideration 
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FIG. 1. (color online). Model schematic, (a) We simulate radial dis¬ 
placement, tracking the fluid-fluid interface (black line) and fluid pres¬ 
sures (increasing from blue to red), (b) Zoom in showing the lattice 
of particles and pores. The interface is represented as a sequence of 
circular menisci, touching particles at contact angle 0 , with curvature 
set by the local capillary pressures. Menisci can be destabilized by: 
(c) burst; (d) touch; or (e) overlap. Brown arrows indicate direction of 
advancement, destabilized arc in dash. 


of details down to the molecular level dim. Considering the 
relative impact of these mechanisms, as well as the complex¬ 
ity and ambiguity involved in their implementation, we em¬ 
phasize viscosity-related mechanisms and exclude liquid films 
and contact line/angle dynamics (see 03 for elaborated dis¬ 
cussion). Consequently, our model provides the coupled ef¬ 
fects of wettability and dynamics in large, disordered domains, 
improving upon existing models which either ignore dynamics 
and/or wettability effects, or are limited by computational cost 
to small domains csiEa. 

We construct a disordered medium by placing cylindrical 
solid particles on a triangular lattice (spacing a), selecting the 
particle diameters d from an assigned distribution; here uni¬ 
form, d e [1 — A, 1 -|- A]J, where A G (0,1) is the degree of 
disorder, and d < a is the mean diameter. The triangular cell 
delimited by a particle triplet defines a pore of volume V, con¬ 
nected to three neighbors via throats of width 2p < a [Fig.[^. 
The fluid-fluid interface is represented by a sequence of cir¬ 
cular arcs (menisci); each arc intersects a pair of particles at 
the prescribed contact angle 0 (iSh , with a radius of curvature 
R ^ 7 /Ap related to the capillary pressure Ap via the Young- 
Laplace law, where 7 is the interfacial tension. 

We consider three types of capillary instabilities ifT^ : (1) 
Haines jump or burst, when the curvature exceeds a threshold; 
(2) touch, when a meniscus intersects a third particle; and (3) 
overlap of adjacent menisci, destabilizing each other [Fig.[ 2 c- 
e)]. Overlap (termed “Melrose event” in CHI) is a nonlocal, co¬ 
operative mechanism affected by the menisci in multiple pores, 
smoothing the interface ifT^ . 

Meniscus instability causes its incipient advancement. Both 
stability and advancement rate depend upon the pressure dif¬ 
ference across each meniscus. Pore pressures and filling rates 


are provided by the fluids’ viscous resistance, evaluated from 
the flow throughout the network of contiguous pores occupied 
by same fluid and through throats with unstable, advancing 
menisci. Flow is resolved via conservation of fluid mass in 
each pore, qj = 0 (summing over all neighboring pores 
j). Assuming Stokes flow, q = CVp provides the interpore 
flow rate, where C ^ p^ j is the conductance and Vp = 
(ft - P)/P- An effective viscosity, peff = {Pi - Pd) ^ + Pd 
allows using q to evaluate both flow of a single fluid between 
two pores and filling rate fH. Here pd and pi are the defending 
and invading fluid viscosities. The filling status of the invaded 
pore, 0 < ^ < 1, is updated according to the inflow from 
throats with unstable menisci q^^^ = qu at each time step 
t, + At) = -f- q^^^{t)At/V. Front readjustments 

are incorporated by considering partially-filled pores, which 
can re-empty upon reversal of meniscus advancement direc¬ 
tion. When pore invasion is completed (^ = 1), the interface 
configuration is updated 03 . The above provides a simple 
description of the invasion dynamics without explicit geomet¬ 
rical calculations of changes in fluid volume from changes in 
menisci curvature. We enforce a constant injection rate from a 
radial region of several pores (inlet), stopping the simulations 
when a boundary (outlet) pore is invaded. 

Our simulations exhibit the experimentally-observed inva¬ 
sion regimes ca: viscous fingering in rapid injection irrespec¬ 
tive of the wettability, capillary fingering in slow drainage, and 
stable, compact displacement in slow imbibition [Fig. [^a)]. 
For a fixed fluid pair (constant Pd/j), the dimensionless flow 
rate is provided by the capillary number, Ca = Pd^l^f, com¬ 
puted from the velocity v = {V tot I Uot) I Aouu where V tot is 
the volume drained during the simulation time Rot, through 
the outlet cross-sectional area Aout- To simulate the injection 
of air into water-glycerol saturated beads in ifTSl . we used the 
following parameters: 7 = 67-10“^ N/m, pi = 1.8-10“^ 
Pa-s, Pd = 5-10“^ Pa-s, a = 500 pm, d = 0.54a, system size 
L = 260a (260x300 particles), and A = 0.81 (providing the 
wide variation of aperture sizes in random bead packs O). 

We characterize the patterns quantitatively via the length of 
the fluid-fluid interfaces (including trapped regions) normal¬ 
ized by the invaded area, I/inter 123 , and the mean finger width 
w (in lattice units a ll28]| ). Viscous fingering is characterized by 
thin fingers of a single pore width, w ^ 1, and long, highly ir¬ 
regular interfaces, I/inter ~ 1, whereas compact displacement 
provides a smooth, rounded front, I/inter 1, with a diverging 
finger width, re ^ 1. In capillary fingering, trapping provides 
long, fractal interfaces, which a patchy, thick pattern composed 
of multiple thinner, contiguous fingers [Fig. |^b-c)]. The ro¬ 
bustness of our characterization is demonstrated by the consis¬ 
tency among four realizations (same particle size distribution). 
Our simulations capture the decrease in finger width with imbi¬ 
bition rate, providing w ^ Ca“^ with z/ « 0.6 [^ = 120°, see 
inset of Fig.j^c)]. While the small difference from u = 0.51 
in CD can be explained by the use of different fluids (and 
0), we note that saturation of re ^ 1 at high Ca exacerbates 
the quality of fit. We also find that sweep efficiency decreased 
sharply between compact displacement, capillary and viscous 
fingering, however non-monotonically C3 
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FIG. 2. (color online), (a) Simulated invasion patterns, characterized 
by: (b) interface length, Linter, and (c) finger width, w. Rapid injec¬ 
tion (high Ca) leads to viscous fingering (VF) with irregular interfaces 
(Winter ~ 1) and thin fingers (w ^ 1). As Ca is decreased, the pat¬ 
terns transition towards capillary fingering (CF) with multiple trapped 
clusters and relatively long interfaces in drainage (low 0), or compact 
displacement (CD, Linter ^ 1, ic 1) in imbibition (high 0). The 
inset of (c) shows a fit of ic ~ Ca“^ for 0 = 120° providing h> ^ 0.6. 
Error bars show the standard deviation among four realizations. 


The crossover between the invasion regimes depends on the 
interplay between three mechanisms: (i) continuous growth of 
thin fingers; (ii) intermittent interface advancement at different 
locations, trapping the defending fluid behind; and (iii) simul¬ 
taneous advancement of large parts of the interface, keeping 
it smooth. Finger growth in (i) is driven by destabilization of 
the entire interface at high Ca, where high defending fluid pres¬ 
sure in the “gulfs” between Angers allows only the Anger tips to 
advance. This screening effect [demonstrated by the pressure 
halo in Fig. [Ja)], where Laplacian-driven growth dominates 
over heterogeneity |[22l [23l, promotes viscous Angering (see 
Videos la and lb in ini)- In (ii), disorder in entry pressures 
leads to capillary Angering at low Ca and 0 (Video 2 in (TTl ). 
At low Ca and high 0, the dominance of overlaps (Fig. en¬ 
hances mechanism (iii), where invasion in one location desta¬ 
bilizes the interface in adjacent pores (Video 3 in (ni, experi¬ 
mentally observed in cax resulting in compact displacement. 
We emphasize that although bursts are the dominant invasion 
mechanism at low 0 irrespective of Ca (Fig. |^, the change in 
driving mechanism (from i to ii as Ca is decreased) leads to 
different patterns. 

We rationalize the invasion behavior by evaluating the mag¬ 
nitude of the forces driving mechanisms (i)-(iii). We predict 
the transition between viscous Angering and capillary flnger- 
ing/compact displacement through a capillary number mod- 
ifled to account for the contact angle, A^ca = ^P±/^P\\- 
Here 6p_\_ is the pressure drop driving growth of individual 
Angers perpendicular to the interface (along the direction of 



FIG. 3. (color online). Occurrence of nonlocal, cooperative pore fill¬ 
ing (number of overlaps out of all instability events) from 208 simula¬ 
tions (gray dots). Increasing 0 enhances overlaps, manifested macro- 
scopically by a more stable displacement. Dashed lines mark the the¬ 
oretical phase boundaries predicted by scaling (see text). 


the externally-applied pressure drop), and 6p\\ is the capil¬ 
lary pressure promoting lateral growth of the interface. We 
evaluate Sp_\_ from the pressure drop in the viscous defending 
fluid over a characteristic length L_\_, 6p±_ ^ Vp^L^, where 
\/p± ^ with permeability k ^ o? and L^ ^ a. We use 

the critical burst curvature Rc (Eq. S2 in (TTl) to evaluate the 
capillary force, 6p\\ ^ ^/Rc, providing 

^Ca = Ca I — P sm‘^0 — IcosO^ , (1) 

where I = d/a is the dimensionless microscopic characteristic 
length. 

For slow injection, we explain the transition between capil¬ 
lary Angering and compact displacement via the “Cooperative 
number” A^coop, a dimensionless parameter evaluating the like¬ 
lihood for pore Ailing by overlaps, 

A^coop = cos ^ — [sin^6> + (2) 

where A^coop = 0 is the geometrical condition for two arcs to 
overlap exactly at their threshold (burst) curvature, such that 
^coop > 0 implies overlap preceding burst ifTTll . Here 0 is the 
local front shape, deflned by the angle between two adjacent 
menisci [Fig. [Je)]. Since the macroscopic pattern is a conse¬ 
quence of numerous invasion events occurring at (p which vary 
in time and space, A^coop represents the instability statistics of 
the entire sample and simulation time: a larger A^coop value im¬ 
plies a higher fraction of overlaps; said differently, for a given 
Ncoop value not all pores will be invaded by the same instabil- 
ity (Fig.[^. Here we compute A^coop using p = 120°, which 
we found to be most representative for our system ifTTl . 

Our scaling analysis predicts the mode of invasion. For rapid 
injection, Nca ^ ^Ca* implies dominance of viscous forces 
leading to viscous Angering. Here, the critical value scales as 
^Ca* ^ ~ 4-10“^, suggesting a dependence on 

the macroscopic characteristic length—the system size E1[71. 
For slow injection, Nca ^ ^Ca*, capillary forces govern and 
invasion becomes strongly dependent on the wettability: for 
nonwetting invasion, A^coop < 0 predicts capillary Angering 
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FIG. 4. (color online). Phase diagrams of immiscible displacement: 
(a) interface length, Linter and (b) finger width, w. At high flow rates, 
Nca ^ predicts viscous fingering (VF) with long, fractal inter¬ 
faces and thin fingers (Linter~1, w^l). At low rates, Nca A^Ca*, 
invasion is controlled by the wettability: for drainage, A'coop < 0 im¬ 
plies capillary fingering (CF), whereas for imbibition A'coop > 0 indi¬ 
cates compact displacement (CD, Linter^ 1, due to smoothing 

by cooperative pore filling. Dashes show phase boundaries from scal¬ 
ing analysis, Nca = A^Ca* ~ 4-10“^ and A^coop = 0. Dots mark 
data from 208 simulations at various Ca and 0. 


viscous instabilities become dominant. Our analysis quanti¬ 
fies the competition between mechanisms governing the dis¬ 
placement stability, insight that could be exploited in technolo¬ 
gies such as microfluidics, hydraulic fracturing and oil recov¬ 
ery lEiiiiiisa. Furthermore, our approach—a set of local rules 
providing a minimal description of the microscopic physics in a 
sufflciently-large domain to capture the emergent macroscopic 
behavior—could provide a new modeling paradigm for other 
problems of front propagation in disordered media, in which 
competition between local disorder, short-range cooperativity 
and global screening play a role, such as active media and spin 
glasses El. 
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caused by disorder in capillary (burst) thresholds, whereas for 
wetting invasion A^coop > 0 implies cooperative motion of 
large parts of the interface (overlaps) and a compact pattern, in 
agreement with our simulations (Fig.|^ and experiments ca. 

The displacement depends on the underlying medium ge¬ 
ometry, including disorder, mean particle size and porosity, in 
a nontrivial manner, as it affects, together with wettability and 
flow rate, the portion of the pore space sampled by invasion. 
For example, more pores would be invaded as disorder and Ca 
are decreased, while decreasing Ca and increasing 0 restricts 
invasion to smaller pores. For the current geometry with high 
porosity (^0.67) and disorder (A = 0.81), A^ca is relatively in¬ 
sensitive to 6, increasing by a factor of ^2.5 from 5° to 120°. 
According to Eqs. the sensitivity of A^ca to 0 increases 

and the threshold angle (corresponding to A^coop = 0, here 
0 = 87°) decreases with particle size 1. 

In this Letter, we have studied the unstable case of high 
disorder and viscosity ratio. Noteworthy perspectives, which 
we intend to study with our model, include the impact of the 
medium geometric properties, viscosity ratio, gravity, and ma¬ 
trix deformations. Our preliminary simulations suggest that de¬ 
creasing the disorder stabilizes the displacement, in agreement 
with Ena. Increased stability is also expected by decreasing 
the viscosity ratio EiEiia or introducing gravity la. Particu¬ 
larly interesting is the coupling with fracturing and particle re¬ 
arrangements, which signiflcantly affects nonwetting invasion 
into granular media fTI I^I^ . 

In conclusion, we elucidate the combined impact of wetta¬ 
bility and dynamics on immiscible displacement in disordered 
media. Our novel model provides the spatiotemporal nonlo¬ 
cal effects of interface dynamics, which are crucial even for 
slow flows due to the intrinsic timescale of interfacial jumps 
which can be orders of magnitude smaller than of the bulk 
flow 1201 . thereby explaining classical yet unresolved observa¬ 
tions. We show that increasing the wettability of the invading 
fluid promotes cooperative pore Ailing that stabilizes the inva¬ 
sion, and that this effect weakens as flow rate increases and 
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